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Abstract 

We point out an equivalence between the discrete velocity method of solving the 
O ' Boltzmann equation, of which the lattice Boltzmann equation method is a special 

o : 

^ ' example, and the approximations to the Boltzmann equation by a Hermite polyno- 

Q^ • mial expansion. Discretizing the Boltzmann equation with a BGK collision term 

c/3 ■ 

d ' at the velocities that correspond to the nodes of a Hermite quadrature is shown to 

p ' be equivalent to truncating the Hermite expansion of the distribution function to the 



o 



c^ 



corresponding order. The truncated part of the distribution has no contribution to 
the moments of low orders and is negligible at small Mach numbers. Higher order 



/\ • approximations to the Boltzmann equation can be achieved by using more velocities 



in the quadrature. 
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The Boltzmann equation is a well accepted mathematical model of a fluid at the microscopic 
level. It describes the evolution of the single particle distribution function, /(x, ^, t), in the phase 
space (x, ^), where x and ^ are the position and velocity vectors respectively. This description 
of a fluid is more fundamental than the Navier-Stokes (NS) equations. It has a broader range of 
application and provides more detailed microscopic information which is critical for the modeling 
of the underlying physics behind complex fluid behavior. However, direct solution of the full 
Boltzmann equation is a formidable task due to the high dimensions of the distribution and the 
complexity in the collision integral. Among the various techniques developed [Jl]], the discrete 
velocity method was introduced [0 based on the intuitive assumption that the gas particles can 
be restricted to have only a small number of velocities. The lattice Boltzmann equation (LBE) 
method formally falls into this category. 

The development of the LBE method for simulation of fluid dynamics was independent of the 
continuum Boltzmann equation. The discrete LBE was first written to describe the dynamics of the 
distribution function in the lattice gas automaton (LGA) [0,0], in which the fluid physics is simu- 
lated at the microscopic level by "Boolean" particles moving with discrete velocities on a regular 
lattice, mimicking the motion of the constituent particles of a fluid. A Bhatnagar-Gross-Krook 
(BGK) collision model [|5|] was later adopted in the LBE in place of the complicated collision 
term [H,^. In this lattice Boltzmann BGK model, the equilibrium distribution is chosen a posteri- 
ori by matching the coefficients in a small velocity (Mach number) expansion so that the correct 
hydrodynamic equations can be derived using the Chapman-Enskog method. 

Recently it has been argued [j^J^] that the LBE method can be derived from the continuum 
Boltzmann equation with a BGK collision model. In the new derivations, the Maxwellian dis- 
tribution is Taylor expanded to second order in the fluid velocity scaled with the sound speed. 
Abe ^ employed a special functional form for the distribution function so that the macroscopic 
fluid variables are completely determined by the values of the distribution function at a set of 
discrete velocities. By noticing that in the Chapman-Enskog calculation, the functional form of 
the equilibrium distribution function in velocity space is only relevant in the calculation of the 
low-order moments, and for the Taylor expanded Maxwellian, those moments can be calculated 



exactly using a Gaussian quadrature, it is concluded that the NS equations can be derived from the 
Boltzmann equation evaluated on the nodes of the quadrature [^. On substituting the weights of 
the corresponding quadrature into the expansion of the Maxwellian, the coefficients of the LBE 
equilibrium distribution function are recovered. The Boltzmann equation evaluated at the discrete 
velocities can then be further discretized in x and t in various ways for numerical integration ^. 
The LBE models are shown to correspond to solving the discrete Boltzmann equations with a 
particular finite difference scheme [[T0|]. 

The recovery of the NS equations from the Boltzmann equation by using a small number of 
collocation points in velocity space is not accidental. Almost half century ago, Grad [ [TT| ] in- 
troduced a sequence of approximations to the Boltzmann equation by expanding the distribution 
function in terms of Hermite polynomials in velocity space. The Hermite coefficients are directly 
related to the macroscopic fluid variables such as density, velocity, internal energy, stress and so 
on. By keeping Hermite polynomials of up to third order, Grad obtained a system of equations 
for thirteen moments of the distribution function. This system of equations, known as the "13 
moment" approximation, was argued to be a better approximation than the Chapman-Enskog cal- 
culation [|T^|T3]] . By noticing that the Hermite coefficients for a given function can be estimated 



using a Hermite quadrature formula, and that this estimation is exact when the function satisfies 
certain conditions, an important correspondence between the LBE method and the approximation 
by Hermite polynomial expansion can be immediately identified. In this Letter, we show that 
by discretizing the Boltzmann-BGK equation at a set of velocities that correspond to the nodes 
of a Gauss-Hermite quadrature in velocity space, we effectively project and solve the Boltzmann 
equation in a subspace spanned by the leading Hermite polynomials. The truncated part of the dis- 
tribution has no contribution to the low-order moments that appear explicitly in the conservation 
equations. 

We start from the following Boltzmann-BGK equation: 

§|- + ^v/ = -i(/-/(°)), (1) 

where, r is a relaxation time, f^^^ is the Maxwellian 
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where D is the dimension of the space, ks is the Boltzmann constant, and m is the mass of the 
molecule. The mass density, p, macroscopic fluid velocity, u, and the temperature, T, are all 
functions of x and t. We introduce the dimensionless quantity 9 = Trrio/Tom, where Tq is a 
characteristic temperature and mo is an unit of the molecular mass. After rescaling the velocities. 



^ and u, in units of the constant cq = w/csTo/mo, which is the sound speed in a gas consisting of 
molecules of mass uiq and at temperature Tq, the Maxwellian takes the following simple form: 

f(o) = P e-^l^-"l' (3) 

The introduction of mo is for the gas mixtures of components with different molecular masses. 
For a single component system we can chose m = mo and have 9 = T /Tq. If the time and length 
scales to and L are chosen so that L/to = cq, the dimensionless Boltzmann-BGK equation will 
have the same form as Eq. ([l|) with r being the dimensionless relaxation time. The mass density, 
p, the dimensionless fluid velocity, u, and the dimensionless internal energy density, e = DO /2, 
are expressed as the velocity moments of the form / f(p{^)d^, with ip = I, ^, and ^'^ respectively: 



P 



j fdi pu = Jf^d$, 2pe + pu^ = J fedt (4) 



In discussions hereafter, the dimensionless variables are used unless otherwise specified. 

Grad introduced the approximations by Hermite expansion in his "13 moment" system [jllp. 
Defining the following weight function 

where ^^ = ^ ■ |, as argued by Grad, if c<j"^/^/ is square integrable, i.e. if / approaches zero faster 
than e~^ /"^ as ^ -^ oo, the following Hermite expansion is valid in the sense of mean convergence: 

OO 1 

/(x,^,t) = ui^) Y: -,at\^.t)'Ht\i). (6) 

where H^'^^ is the nth order Hermite polynomial [[1^]. This expansion, also known as the Gram- 
Charlier series, was also used in the solution of the Vlasov equation [[TS]]. Clearly, the Maxwellian 



in Eq. (|]) has such an expansion if the choice of Tq and niQ ensures that 9 < 2. On the right- 
hand-side of Eq. (^, both a^"^ and H^"'^ are symmetric tensors of order n. The subscript i is an 
abbreviation for the multiple indices {ii,- ■ ■ , in}, and the products denote contraction on all the 
n indices. The Hermite polynomials are a set of complete orthonormal basis of the Hilbert space 
with the inner product (/, g) = J Lufgd$,. They satisfy the orthonormal relation: 

Ju;nt^n^;^d^ = 6,nnS.j, (7) 

where 6ij = 1 if i = {ii, ■ ■ ■ , im} is a permutation of j = {ji, ■ ■ ■ ,jn}, and Sij = otherwise. 
For any function /, the nth Hermite coefficient can be obtained by the following equation: 

/oo „(m,) „ 
fn^''\^)d^ = E ^ / ^'n'r^'nfh^ = a^^\ (8) 

where the factor m! is the number of terms in 8ij. The function / is completely determined by all 
of its Hermite coefficients. 

The moments given in Eqs. (Q) are invariants of both the original Boltzmann collision term and 
the BGK collision model. The hydrodynamic equations are simply the corresponding conservation 
equations: 



d_ 
di 



J f^d^ + V- J f^^d^ = 0. (9) 



The lowest order moments have the most significant contribution to the macroscopic hydrodynam- 
ics. Since the Hermite expansion has the feature that a velocity moment of a given order is solely 
determined by the Hermite coefficients up to that order and are not changed by the truncation 
of the higher-order terms, a sequence of approximations to Eq. (jl]) can be made by seeking the 
approximate solution of the following form: 

^ 1 
/(x,^,t) = a;(0 E -a^\^,t)n^\^). (10) 

n=0 '^• 

The momentum and energy conservation equations explicitly involve moments of up to the second 
and the third order respectively. It is necessary to require A^ > 2 if the momentum equation is to 
be obtained and A^ > 3 if the energy conservation equation is needed. 



By the approximation above, we have assumed that the distribution function lies entirely in the 
subspace spanned by Hermite polynomials up to order A^. For higher orders, 

a(") =0 if n> N. (11) 

Although the terms that are truncated do not appear explicitly in the conservation equations, they 
affect the fluid variables through their contributions to the dynamic equations of the lower order 
moments. We will return to the validity of the assumption Eq. ( [TT] ) later. 

Let ^j and Wj, z = 1, ■ ■ ■ , d, be the nodes and weights of a quadrature of degree 2N, i.e., if p(^) 
is a polynomial with a degree not greater than 2N, we have 

fuJ{^)p{Od^ = J:^^p{^^)■ (12) 

4=1 

Because fTi^'^'^ /uj is such a polynomial if n < A^, the Hermite coefficients of / can be calculated 
using the values of / at the nodes ^^ as the following: 



~f^in),, ^nj,m(-\^. 






(13) 



where, /, = /(x, ^j,t). The knowledge of /j as functions of position and time is equivalent to 
that of the truncated distribution function itself and therefore, that of the fluid variables calculated 
from the truncated distribution. These variables are: 

^ Wifi ^WiMi 



2pe + fm^ = j: ^. (14) 

By defining the auxiliary variables gi = Wifi/uj{^^), Eqs. ( |T4| ) can be put into a more efficient 
form for computation: 

d d d 

P = ^9i, Pu = ^ gi$,i, 2pe + pu^ = J2 9i^i- (15) 

j=l «=1 4 = 1 

This has the same form as how the fluid variables are calculated in the LBE method, where the dis- 
tribution function is defined from the beginning as the populations of particles moving at discrete 
velocities. 



We now turn to the discussion of the equations that the functions fi satisfy. By directly evalu- 
ating Eq. ([1|) at the nodes ^j, we have: 

^+^.-v/. = -i[/.-/(°)(ej]. (16) 

Because f^^^ has non-zero Hermite coefficients at all orders, on substituting /'•"-'(Cj) into the right 
hand side of Eqs. (p^, the equalities hold only approximately. Clearly the mass density, momen- 
tum and total energy will not be exactly conserved by the collision term on the right hand side 
of Eq. (jT^. For the conservation laws to hold exactly, f^^^ has to be projected into the subspace 
in which / lies. Namely, in Eq. ([T^, the Maxwellian has to be replaced by its A^th order Her- 
mite expansion. Denoting the values of the auxiliary variables gi corresponding to the truncated 
Maxwellian by gl , Eq. ( [T6[ ) can be written as 

^ + ^rV^. = -^((7.-#)), t = l,---,d. (17) 

It is to be noted that the positivity of the distribution function is lost in this truncation. 

The first a few Hermite coefficients of the Maxwellian can be obtained using Eq. (Q). They are: 

a(°) = p,ai^^ = pui,afj' = pUiU j + {6 -1) p6ij, and a[fi^ = pUiUjUk+{6-l)p{ui6jk+Uj6ik+UkSij). 
Using the explicit form of the Hermite polynomials, at the second and the third orders, we have 



~(2) 

9i = Wip 



~(3) 

91 = WiP 



1 u^ 

l + 7. + eru + -(^ru)'-y 

1 u^ 



(18) 



+ \{i^■^f-\{i.■^)u' 



(19) 



where 7^ = {^j - D){e ~ l)/2, Q = {^f -D-2){e- l)/2. Eqs. (|T7|)-(|19|) are the projection of 
the Boltzmann-BGK equation in the subspace spanned by the leading Hermite polynomials. They 
are in the configuration space (t, x) and have a linear differential operator on the left-hand side. 
The fluid variables defined by Eqs. ( [T3| ) obeys the NS hydrodynamics. As previously shown ^]^, 
the LBE equilibrium distributions of Refs. [§,0] are obtained when the proper nodes and weights 
are substituted into Eq. ([T^), and the lattice Boltzmann equations are particular finite difference 
discretizations of Eq. ([T7|). 



It can be easily verified that moments of up to second and third orders calculated from Eqs. (jT^) 
and (|T9|) respectively are those of the Maxwellian. In particular, the tensor, p^ = / f^^^^i^jd^ = 
pUiUj + p95ij, survives the truncation. The hydrostatic pressure is given by the dimensionless 
equation of state, p = p9, which translates to p = -^ksT in laboratory units, and the sound speed 
is \/9 as expected. When measured in the magnitude of one of the nodes of the quadrature, e.g. 
^i, as in the LBE models, the sound speed is VO/^i. In a single component isothermal system, 9 
becomes a free parameter which can be used to adjust the nominal sound speed with respect to the 
nodes of the quadrature. When 9 = 1, the Maxwellian has a very simple expansion: 

/(o)=e^py i^u(")7^("). (20) 

and the truncated part of the distribution function is proportional to the power of the Mach number. 
Eqs. (p^-(P^ are also simplified since 7j = C* = 0. For a multiple component system, a different 
9 ~ 1/m has to be chosen for each component if all the components are at thermal equilibrium. 
This requirement was found necessary to obtain a correct equation of state in a previous multiple 
component LBE model [p^. 

It should be noted that the truncation made in Eq. ([n|) is similar to, but not exactly the same 
as the third order approximation in the Grad "13 moment" system. In the latter the distribution 
function is expanded in velocity space around the local fluid velocity before it is truncated. This 
is certainly a better approximation than expanding around the absolute equilibrium, and corre- 
sponding computational methods can be developed. However, it is not possible to use such an 
expansion in the method discussed above because it would yield a set of nodes that depend on the 
local velocity. 

The difference between the two expansions can be estimated by expanding the following ap- 
proximated distribution function of the "13 moment" system [ ]TT| ] around the absolute equilibrium: 

/(°)E^^«(^-u). (21) 

i=o '■■ 

where, 6'^°) = 1, fe^^^ = 0, and bl^ = 0. The Hermite coefficients in the expansion around the 
absolute equilibrium, a^"\ can be calculated as the following: 



a(") 



Jfio) J2 ^7^W(^ - u)n^-\^)d^ (22) 

i=0 

pE ^ f^im^Hm^'^K^ + u)rf^ (23) 



Noticing that 7^(")(^ + u) = Er=o u("~^^H(^)(^), we find that a(°) = p, a^^) = pu, a^^) = 
p(u(2)+6(2)),andforn>3, 

With non-zero Hermite coefficients at all orders, the distribution function in the "13 moment" 
system does not meet the assumption in Eq. dTT]). Since these coefficients are proportional to the 
power of the Mach number, Eqs. (|19|) approximate the "13 moment" system only at the small 
Mach number limit. This approximation can be improved by using a quadrature formula with a 
higher degree. 

Except for one dimension, finding the quadrature formula with minimum number of nodes 
for given geometry, weight function, and degree of accuracy is generally an unsolved problem. 
As the minimum requirement, a formula of 4th degree is needed for isothermal systems and 6th 
degree for thermodynamic systems. For the weight function in Eq. (^, quadrature formulae of 
different degrees are listed in Ref. [[T7|]. Some of those are believed to be minimum without proof. 
In two dimensions, the minimum formula seems to be the 4th degree, 6-point formula (origin 
and the vertices of a pentagon) for isothermal models and the 7th degree, 12-point formula for 
thermodynamics models. In three dimensions, the minimum formulae are those of 5th degree, 13- 
point (origin and the vertices of a regular icosahedron) and 7th degree, 27-point for isothermal and 
thermodynamic systems respectively. The nodes of these formulae usually do not coincide with 
that of a regular lattice. Eqs. ( [T7| ) have to be solved using schemes such as the finite difference 
method \M^. 



The methodology of the discrete Boltzmann equation can be summarized as the following. 
The discretization of the continuum distribution function into values at the nodes of a quadrature 
formula is equivalent to the truncation of the high-order terms in the Hermite spectral space. The 
information that is lost in this procedure is represented by the high-order Hermite polynomials 



which do not explicitly appear in the conservation equations. This error is negligible at small 
Mach numbers and can always be made smaller by using a quadrature of a higher degree. With 
such a discretization, the Boltzmann equation becomes a homogeneous set of linear equations in 
the configuration space. Comparing with the non-linear NS equations, these equations are easier 
to solve, have a broader range of application, and more importantly, allow the underlying fluid 
physics to be simulated directly at the cost of a macroscopic simulation. In addition, higher order 
approximations to the Boltzmann equation can be easily achieved by adding more points to the 
system. 

Some of the limitations that LBE methods inherited from the Boolean LGA models can be 
removed with the present formulation. The equilibrium distribution is now obtained through a 
systematic orthogonal expansion of the Maxwellian, eliminating the tedious parameter-matching 
procedure which usually produces results that are not unique and contain erroneous terms at higher 
orders. The inflexible lattice structure and time stepping scheme of the LBE method are inconve- 
nient for practical applications and often result in poor stability. By realizing that the LBE models 
are merely simple and rather primitive finite difference representations of the discrete Boltzmann- 
BGK equation, we can employ more sophisticated numerical techniques in solving these equations 
with better efficiency, stability and flexibility. We defer the detailed discussion and numerical ex- 
amples to a future publication. 

The authors thank Dr. Gary Doolen and Dr. Nicos Martys for helpful discussions. 
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